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For inhomogeneous systems with interfaces, the inclusion of long-range dispersion 
interactions is necessary to achieve consistency between molecular simulation calcu- 
lations and experimental results. For accurate and efficient incorporation of these 
contributions, we have implemented a particle-particle particle-mesh (PPPM) Ewald 
solver for dispersion (r~ 6 ) interactions into the LAMMPS molecular dynamics pack- 
age. We demonstrate that the solver's 0(N log N) scaling behavior allows its appli- 
cation to large-scale simulations. We carefully determine a set of parameters for the 
solver that provides accurate results and efficient computation. We perform a series 
of simulations with Lennard- Jones particles, SPC/E water, and hexane to show that 
with our choice of parameters the dependence of physical results on the chosen cutoff 
radius is removed. Physical results and computation time of these simulations are 
compared to results obtained using either a plain cutoff or a traditional Ewald sum 
for dispersion. 
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I. INTRODUCTION 



Despite their weak r -6 scaling, dispersion forces "play a role in a host of important phe- 
nomena such as adhesion; surface tension; physical adsorption; wetting; the properties of 
gases, liquids, and thin films; the strength of solids; the flocculation of particles in liquids; 
and the structures of condensed macromolecules such as proteins and polymers. Unsur- 
prisingly, their contributions to intermolecular interactions are accounted for in the vast ma- 
jority of the nonbonded potentials applied in molecular simulations. Typically, dispersion 
interactions are only considered within a cutoff of around 1 nm. For homogeneous systems, 
the contributions of long-ranged interactions can be estimated efficiently and accuratelyP 
For inhomogeneous systems, however, these corrections are inaccurate, and computational 
requirements has precluded the inclusion of long-range dispersion interactions, even though, 
as can be seen from the applications above, they are especially relevant for these systems. 
The necessity of incorporating the long-range effects of dispersion forces has been shown in 
numerous studies on surface simulations^^ of which only some are referenced here, but also 
in simulations near the critical pointp^and in simulations of protein-ligand binding.^ 

Various correction methods have been proposed for incorporating long-range dispersion 
contributions. Most molecular simulation packages already include energy and pressure 
corrections assuming homogeneous systems. Similar "on-line" methods that can be applied 
during simulations have been presented by Guo et al.^for Monte Carlo and by Mecke et alP 
and JanecekPfor molecular dynamics (MD) simulations. Chapela et alPand Blokhuis et alP 
have developed a tail correction for simulated surface tensions that can be added after the end 
of the simulations. The aformentioned correction methods are applicable only to simulations 
with planar interfaces. The use of a twin-range cutoff has been proposed by Lagiie et alP^l 
Wu and Brooks^ present the isotropic periodic sum for electrostatic interactions, but the 
method can be extended to incorporate long-range dispersion forces. 

An alternative to these "correction" -based schemes is to include the long-range interac- 
tions explicitly using Ewald summation, which was originally developed for the treatment 
of Coulomb forces.^ This method was developed for dispersion by Williams,^ Perramp^ 
and Karasawa and Goddardp*' and later applied to surface simulations by Lopez-Lemus et 
al.P Grest and co-workers JEHUl Ou-Yang et al.™ and Alejandre and ChapelaP^ Among the 
previously mentioned methods for treating long-range dispersion interactions, the Ewald 
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sum is usually the most accurate, most reliable, and most versatile. However, its 0(N 3 ^ 2 ) 
scaling' prohibits its application in large-scale simulations. This problem can be overcome 
by applying grid-based Ewald summation methods,^^ w hose scaling, because of the use 
of fast Fourier transforms (FFT), is 0(N log N). While frequently used for Coulomb inter- 
actions, such methods have also been applied to dispersion interaction by Essman et alP^ 
and Shi et alP^ In the former work, the dispersion part is only adressed very briefly for 
particle- mesh Ewald (PME), and not for PPPM. The latter puts a stronger focus on the 
PPPM and dispersion interactions. We feel, however, that the description is incomplete; 
for example, the equations for the energy and virial and the exact formulation of the true 
reference force are not given. Furthermore, we provide a simpler method for calculating 
the pressure tensor required for calculating surface tensions, and outline reasons why their 
simulated surface tensions do not agree with other reported results for SPC/E water. 

We present the results of the development and implementation of a particle-particle 
particle-mesh (PPPM) solver for r~ 6 dispersion interaction in the LAMMPSP^ molecular 
dynamics package. The theory is given in Section [TT} Error estimates for the real and re- 
ciprocal space contributions, as well as a discussion on the limits of the error estimate, are 
presented in section |III| The scaling behavior of the developed algorithm is presented in 



Section IV We have performed a variety of interfacial simulations, as long-range dispersion 
interactions are known to have a significant effect on simulated surface tensions. The theory 
for calculating surface tensions is briefly reviewed in Section |Vj A set of parameters for 
performing successful simulations with the PPPM for dispersion is determined in Section 



VI. We use these parameters in Section VII to study the surface tension of Lennard- Jones 



(LJ) particles, hexane, and SPC/E water. Section VIII contains a brief comparison of the 
simulation time of different solvers. We summarize our findings in the final section. 



II. FORMULATION OF THE MESH-BASED DISPERSION SUM 

Excellent reviews on traditional and mesh-based Ewald sums are given by Hockney and 
Eastwood,^ Essmann et al.,^ and Deserno and HolmP^ Karasawa et al.^ provide a com- 
prehensive description of the traditional Ewald sum for dispersion interactions. To make 
the presentation as self-contained as possible, we provide a complete summary of the PPPM 
algorithm as applied to dispersion interactions, compiled from the above references. 
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Particle Distance 

FIG. 1. Total potential and split potential for r -6 . When not using the Ewald sum, the whole 
area under the blue curve is solved in real space, whereas only the area under the green curve is 
solved in real space when using the Ewald sum. The error in the calculation is related to the area 
under the curves beyond the cutoff. Using the Ewald technique is thus more accurate. 

Because of its physical origin in the overlap of electron hulls, the repulsive (often r~ n , 
where typically n > 9) part of pair potentials is very short-ranged and can be neglected 
beyond a typical cutoff length of 1 nm. We therefore exclude the repulsive term from further 
consideration. The attractive part between two atomic sites of some commonly used pair 
potentials, such as the LJ or Buckingham potential,^ can be expressed as 

C ■ 

Mattr(^i) = -— (1) 

where r is the distance between particles % and j and Cy is the dispersion coefficient de- 
scribing the strength of the interaction. The goal of the Ewald summation is to split this 
potential into a fast-decaying potential, whose contribution can be neglected beyond a cut- 
off, and a slowly decaying potential, whose contribution can be accounted for in Fourier 
space, as shown in Figure [TJ Its calculation requires a Fourier transform of the dispersion 
coefficients into the reciprocal space. 

The benefit of mesh-based Ewald methods, such as PPPM, is that the dispersion coeffi- 
cients are distributed onto a grid, which permits application of the FFT for the calculation 
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of the dispersion coefficient density. The calculation of the mesh Ewald sums in reciprocal 
space requires additional steps. The dispersion coefficients have to be distributed onto a 
grid and transformed into reciprocal space to calculate their interactions in Fourier space. 
The resulting potential must then be derived and backinterpolated onto the atomic sites to 
obtain the forces. 

The dispersion energy of a system with the potential above is given by™ 

keM 

id * 

where /3 is the Ewald parameter for dispersion interactions, is the distance between 
particle % and the nearest image of particle j, and V is the volume of the simulation box. G 
is the optimum dispersion influence function, which has a different form than the electrostatic 
influence function. <S| is a function of the location and strength of the dispersion sites. The 
k vectors form the discrete set {27rn/L}, where L is the length of the box vectors and the 
components of n = (n x , n y ,n z ) are integer values that are zero for the center node of the 
grid. The first sum in equation [2] is over all pairs of atoms. However, as the potential decays 
quickly with increasing interparticle distance, it is only evaluated for particles whose is 
smaller than a chosen cutoff. The second sum is over the reciprocal vectors of all grid points. 

The expression for the optimum influence function, which minimizes the error in the 
calculated forces, i d 21 ' 24 1 

d = D(k)S^ P &Mk + ¥m)R (k+frn^ (3) 



D(k)l 2 £ me z. U 2 (k + fm) 



with 



UQt) = ^y^, (4) 



where P is the interpolation order and 



w(-P) _ ^3 ( sin (kxh/2) sm(k y h/2) sm(k z h/2) \ 
k x h/2 kyh/2 k z h/2 J 



(5) 



is the Fourier transform of the interpolation function W^ p \ which is described, for example, 
in Ref. [2H and h is the grid spacing. D is the Fourier transform of the differentiation 



operator required for the force calculation. In this study, we use differentiation in Fourier 
space 

D(k) = zk. (6) 

R is the Fourier transform of the true reference force and can, for dispersion interaction, be 
calculated a d 23 ' 24 1 

vr 3 / 2 /? 3 r 9 i 

R(k) = tk-^-f— (1 - 2b 2 ) e~ b + 26 3 v^erfc (6) , (7) 

with b = |k|/2/3. 

For our choice of U, the denominator in equation [3] can be evaluated analytically, as 
shown by Hockney and EastwoocP^or in a more explicit form by Pollock et al.^ The sum in 
the numerator is usually sufficiently converged when terms with |m| < 2 are included. As 
the influence function is independent of the particle coordinates, it needs to be calculated 
only at the beginning of a simulation or when the volume has changed. 

When the dispersion coefficient of a pair of atoms can be expressed using a geometric 
mixing rule, 

Cij = \J CuCjj = QCj, (8) 
as in, for example, the OPLS potential)^ the function <S|(k) can be expressed as 

4 2 (k) = S 6 (k)£*(k), (9) 

where S$ is the complex conjugate of the structure factor S&, which is the discrete Fourier 
transform of the dispersion coefficient density cm on the grid points: 

5 6 (k) = c M (r p ) exp (-zk • r p ) . (10) 

r p GM 

When using an LJ potential, the dispersion coefficients of unlike sites are often determined 
via the Lorentz-Berthelot mixing rule as 

0"i + o"., 



Cy = V^ ^ , (11) 



where e and a are LJ parameters. Equation [sj cannot be used in this case; instead, «S|(k) 
must be calculated by^ 
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4 2 (k) = ^5 6 , fc (k)5* ife (k), (12) 



i=0 



with 

S 6>k (k) = ^ c M/( r p) ex P ■ r p) > ( 13 ) 

r p GM 

where c^m is the dispersion coefficient density on the mesh points obtained from interpo- 
lating the dispersion coefficients 



^-kvGD* (14) 

onto a grid. Because of their symmetry only four of the seven addends in equation 12 have 
to be calculated. Although the imaginary part of each addend is not necessarily zero, the 
imaginary parts of the sum will cancel out identically. 

Additional steps are required for calculating the mesh-based forces. For ik differentiation, 
the dispersion field can be calculated as^l 

E(r p ) = -FFT (ik^G^j (r p ), (15) 

where F~FT indicates the reverse fast Fourier transform. The total force on particle i, based 
on the the real and the reciprocal part, can then be calculated a d 19 * 24 * 



ij ij ij ij 



X 



exp (-r 2 /3 2 ) Tij + CiJ2 E(r p )^( ri - r p ) 



r P er 



(16) 



It should be noted that equation 15 and the second term in equation 16 have to be calculated 
for each of the seven structure factors when the Lorentz-Berthelot mixing rule is used. 
The instantaneous stress is given by^ 



6 6/3 2 3/3 4 f3 6 



ij \ ij ij ij ij- 



exp (-r 2 /3 2 ) t^tw + ^ G(k)S 2 6 (k) 
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3 26 3 v ^Ferfc(6) - 26 2 e ; 
aP ~ fc* 2& 3 v^erfc(&) + (1 - Wfir* a fi 
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where 5 a/ 3 is the Kronecker delta. 



III. FORMULATION OF AN ERROR MEASURE 



Several parameters can be tuned to influence the accuracy of the dispersion PPPM: the 
chosen cutoff radius r c for the sum in real space, the Ewald parameter /3, the grid size, 
and the order of the interpolation function for distributing the dispersion coefficient onto a 
grid. The qualitative influence of the parameters can be understood easily. The real space 
error arises from truncating the pair potential. Increasing the cutoff radius or the Ewald 
parameter, which leads to a faster decaying real space potential, increases the accuracy in 
real space. The precision in reciprocal space depends on the Ewald parameter, the grid 
spacing, and the interpolation order. Decreasing either of the first two or increasing the 
latter of these parameters will lead to higher accuracy in the reciprocal space contribution. 

To choose the tunable parameters effectively, a more quantitative understanding of the 
parameters' influence is required. The following sections present estimates for the error of 
real and reciprocal space contribution to the forces. 

A. Error measure for the real-space contribution 

To describe the real space error, we extend the derivation of Kolafa and PerramP^ for 
Coulomb interaction to r~ e potentials. The sum of the square of the real-space contribution 
of the dispersion interaction of the particles beyond the cutoff r c on a single particle can be 
expressed as 



Assuming that the particles are randomly distributed beyond the cutoff, the sum can be 
replaced by an integral to arrive at 




exp (-2r^ 2 ) r%. 



(18) 




exp (-2r|/3 2 ) 47rdr. 



Using 
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which is valid for B > 0, and 



d(f(x)/2Bx) 
dx 



<f(x) 



for x > A, we arrive at 

* ^2 2 7T/3 10 f 6 6 

AF = c > c — — 1 h 

exp (-2r 2 /3 2 ) , 
which leads to the averaged error in the force 



6 6 
+ 



where iV is the number of particles and 



c?. 



r 2 /3 2 



+ 1 x 



j 

exp (-r 2 /3 2 ) 



+ 1 x 



(19) 



(20) 



B. Formulation of an estimate for the reciprocal space error 

The following sections present an estimate for the error of reciprocal space contribution 
to the forces that is an extension to r -6 potentials of the analytical error measure by Deserno 
and HolmpH for Coulomb interactions. 

Following the reasoning from Deserno and Holmj^ the error in the forces in reciprocal 
space can be expressed by 

^-freciprocal = C y jy^7' (21) 

where Q can, for the optimum choice of the influence function, be calculated as 



R ( k + — m 



keM kmGZ 3 

|D(k) E me za U 2 (k + f m) R* (k + f m) 



(22) 



|D(k)| 2 [En^ 2 (k+fm)_ 

In the following, we will derive an approximation for Q that can be rapidly calculated. This 
approximation is restricted to cubic systems with the same number of mesh points N m in 



each direction and the ik differentiation scheme employed in this study. It is based on the 
assumption that h/3 is small. 

Like Deserno and Holmj^ we exploit the fast-decaying form of R to make the approxi- 
mation that it is sufficient to retain only |m| = in the inner sums containing R. Following 
the same steps, we thus arrive at 



Q 



7T/3 



6 r°° 
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2(P + m) + l dr ' 



(23) 



where cL ^ are coefficients given in Table This equation corresponds to Equation (32) in 
Ref . [3T1 before changing the sum to an integral. Performing the integration leads to the final 
result 

p-i 



Q 



3 

12 



m=0 



h0\ 2{P+m) 2 
2 J 2 (P + m) + 1 



x 



{ Tl + T 2 + T 3 + T 4 + T 5 + T 6 } 



(24) 



with 
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T 2 

T 3 
T 4 
T 5 
T 6 



[2 (P + m) + 3]!! 

" 71 ' 
[2 (P + m) + 5]!! 

" 71 ' 
[2 (P + m) + 7]!! 

471 ' 
2 P+m+3 [2(p + m + 3)]!! x P 4 , 

_ 2 p+ m + 3 [ 2 (p + m + 4)]i! x p 5; 

2-P+m+3 



2 (P + m + 5) + 1 
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where x\\ is the double factorial function 

x!! = x(x-2)!!,(0)!! = (-1)!! = 1, 
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TABLE I. Coefficients required for the calculation of the reciprocal space error estimate. (Reprinted 
from Ref. EH) 
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and V\ is given by 

P+m+l-l 

-P, = 1 _ ,/9 V" 

(2i)\\2 i 



P+m+l-l /„ . 1 n ■■ 



i=0 



C. Numerical Tests 

We performed test runs to examine the accuracy of the real space and reciprocal space 
error estimates. We placed 2000 L J particles randomly in a box with box length 15 a in 
each direction to create a bulk system. In order to test the error estimates for surface 
systems, we placed 4000 LJ particles randomly in a 30 x 30 x 10a 3 box and extended the 
length of the shortest box edge to 30 a afterwards without changing the particle coordinates. 
We calculated the real and reciprocal space forces on the particles for these configurations 
seperately using different values for the Ewald parameter, the grid size, the interpolation 
order, and the real space cutoff. Interpolation orders P = 3, . . . , 6 and 2 k mesh points, 
k = 2, . . . , 6 in each direction were used. Real-space cutoffs of 2.0, 3.0, and 4.0 a were used. 
The error in the forces is calculated as 



AF 



1 N 

_^ (F PPPM_ Fr ct )2; (26) 
i=l 

where Ff PPM is the force calculated with the PPPM and F| xact is the "exact" force calculated 
with an Ewald summational in which we used a large cutoff and a large number of recirpocal 
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FIG. 2. Real space error in the forces measured and estimated with equation 20 for (top) bulk 
system and (bottom) interfacial system. Measured errors are depicted as solid lines, estimated 
errors as dashed lines. From top to bottom, the real space cutoff is (blue) 2. 0<r, (green) 3. Oct, and 
(red) 4. Oct. The estimate works well for the bulk system, but fails for the interfacial system. 

vectors to ensure proper conversion. 

The results of the real space error estimate are given in Figure [2} Results for the reciprocal 
space error estimates are shown in Figure [3j Except for small values of /3, the real space 
error estimate works well for the bulk system. In contrast, the error is underestimated 
in surface simulations. For bulk phase systems, the reciprocal space error estimate with 
equations[2T]and[22]provides very good results. The approximation with equation 24 strongly 



overestimates the reciprocal space error when the assumption that h/3 is small is violated. For 
the interfacial system, the error estimates underpredict the simulation error. Yet, as can be 
seen from these figures and from Figure |6j the error estimates can be useful for determining 
the value of the Ewald parameter for which the accuracies in real and reciprocal space are 
equal, if this information is needed. 

The results shown above demonstrate that the error estimates presented here should only 
be applied to homogeneous bulk systems. Additionally, it should be noted that the error 
estimate for the real space contribution assumes that the errors in the forces partly cancel. 
This cancellation of errors cannot occur for the real space contribution to either energy 
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FIG. 3. Measured and estimated reciprocal space errors, for (top) a bulk system and (bottom) an 
interfacial system. The graphs on the left show increasing interpolation order from P = 3 at the 
top to P = 6 at the bottom, with a fixed grid containing 32 grid points in each direction. The 
graphs on the right show grid density increasing from 2 2 points in each direction at the top to 2 6 
points at the bottom, for fixed interpolation order P = 5. Measured errors are depicted as solid 



lines. Errors estimated with equations [21] and 24 are depicted as dashed lines, those estimated 



with equations 21 and 22 are depicted as dotted lines. 



or pressure. This can be easily seen from the following example: Consider three equal, 
collinear particles, with particle 2 equidistant between particles 1 and 3. The distance 
between particle 2 and the other particles is larger than the chosen real-space cutoff, so that 
none of the real-space forces, energies, or pressures are calculated. If the chosen cutoff radius 
were larger, so that the interactions should be calculated, the force on particle 2 would be 
zero, because the contributions from particles 1 and 3 cancel. The energy and the pressure 
that are exerted on particle 2, however, do not cancel but instead are additive. The reason 
for this behavior is that dispersion interactions, unlike Coulomb interactions, are always 
attractive. Contributions to the energy thus always have the same sign. As the distance 
vectors and force vectors for pairwise interactions always point in opposite directions, the 
contributions to the diagonal components of the virial tensor always have the same sign, too. 
This in turn means that pressures and energies can be underestimated, even when forces are 
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calculated accurately. 

Thus, usage of the above error estimates to set the Ewald and grid parameters is only 
recommended for bulk systems in which neither the energy nor the pressure is relevant, such 
as in simulations for determining diffusivities. We show how to determine parameters for 
interfacial simulations in Section IVIl 



IV. SCALING OF THE ALGORITHM 

The main benefit of mesh-based Ewald methods over traditional Ewald sums is the im- 
proved scaling behavior of the mesh-based approach. To examine the scalability of the 
implemented solver, we have performed simulations with 2™ x 10 3 LJ particles, where 
n = 0,1,... ,10, with the dispersion PPPM solver and the Ewald summation. The den- 
sity was 3.64<t -3 in all simulations. The boxes were always cubic. An energy minimization 
and equilibration over 50 000 timesteps in the NVT ensemble at a reduced temperature 
T* = 0.85 was followed by a simulation over 1000 timesteps in the NVE ensemble. The 
simulation time of the last 1 000 timesteps was used to measure the performance. These 
simulations were executed on a single core of an Intel Harpertown E5454 processor with 
eight 3.0 GHz Xeon cores. 

Automated parameter generation was applied in simulations with the Ewald sumP The 
mesh parameters for the PPPM were set using the error estimate presented in the previous 
section in the following way. The real space cutoff was chosen as 3.0 a. The real space error 
estimate was then used to set the Ewald parameter to obtain a desired accuracy of 0.01 e/a 
in the calculated real space forces. The interpolation order was set to P = 5. Using these 
data, the grid spacing was chosen in a way that the accuracy of the reciprocal space forces 



was smaller then 0.01 e/a by using equation 24 As the conditions for the validity of the 



error estimates are fulfilled for the chosen simulations, the comparison we draw here is for 
different system sizes run with the same accuracy. 

As can be seen from Figure |4| which shows the computation time per timestep, the 
dispersion PPPM approaches the expected scaling behavior of 0(N log N) with increasing 
numbers of particles. Its performance becomes several magnitudes faster than the traditional 
Ewald sum and is thus far more suitable for large-scale simulations. The comparison between 
the different solvers drawn here should be considered qualitative, as we did not examine 
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FIG. 4. Scaling of the Ewald sum and the PPPM method for dispersion interactions 
whether the two different solvers were run with the same accuracy. 



V. DETERMINATION OF SURFACE TENSIONS FROM MD 
SIMULATION 

As the need for incorporating long-range dispersion is especially acute for interfacial 
systems, we have run simulations with explicit interfaces on LJ particles, SPC/E water, and 
hexane to test the efficiency and accuracy of the dispersion PPPM algorithm. This section 
briefly summarizes the method applied to simulate surface tensions. Surface tensions can be 
obtained from MD simulations via two-phase simulations. We use the approach, developed 
by TolmarP^l and Kirkwood and Buff in which the surface tension is expressed via 

i r°° 

7 P =2 J {px{z)-p\\{z))dz, (27) 

where p±(z) = p z (z) is the pressure component normal to the surface and p\\(z) = (p x (z) + 
p y (z))/2 is the pressure component parallel to the surface. Replacing the integral with an 
ensemble average leads to 



i P = y {p±-p\\) = y 



(Px) + (Py) 



(28) 



2 

where L z is the box dimension in the z-direction. The outer factor of 1/2 takes into account 
that the simulated system contains two interfaces. 
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If a cutoff is introduced for the pair potential, the surface tensions calculated with Equa- 
tion [28] will underestimate the correct surface tension of the simulated material. This error 
can be estimated by adding a "tail correction" 7 tai i to the simulated surface tension to 
provide a better estimate of the correct surface tension 

7 ~ lp + Ttaii ( 29 ) 



from the simulation. The correction can be calculated as^ 

x [p{z)p{z — sr) — [pq{z)) ) drdsdz, (30) 

where U(r) is the pair potential, g(r) is the radial distribution function, p(z) is the simulated 
density profile, r c is the cutoff radius for the pair potential, and Pq{z) is the Gibbs dividing 
surface 

p G (z)= p c +^sgn(z), (31) 

where p c is the mean and Ap is the difference of the densities of the coexisting phases. g(r) 
was assumed to be unity beyond the cutoff in the calculations of the tail correction. The 
values for p c and Ap, which were also used to calculate the liquid and vapor densities in this 
study, were obtained from fitting an error function to the simulated density profiled 



VI. INFLUENCE OF THE EWALD AND GRID PARAMETERS ON 
PHYSICAL PROPERTIES 

The parameters used by the dispersion PPPM have a strong influence on both the effi- 
ciency and the accuracy of the simulations. As the presented error estimates fail to describe 
systems with interfaces, we have run test simulations to determine a set of parameters that 
can provide both accurate results and acceptable performance for interfacial simulations. 
These parameters were determined for Lennard- Jones particles and hexane, a nonpolar fluid 
whose intermolecular interactions are dominated by dispersion. Hexane was modeled using 
the OPLS-AA 133 force field. 

Simulations with hexane contained 689 hexane molecules that were placed using PackmoP^ 
in a subvolume around the center of the box with volume 50 x 50 x 150 A 3 . After an energy 
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minimization with a soft potential and several runs with restricted movement of the parti- 
cles, the simulations were equilibrated for 1 000 000 timesteps with a timestep of At = lfs. 
The temperature was set to T = 300 K using a Nose- Hoover^ thermostat with a damping 
factor of 0.1 ps. A PPPA/P with a real space cutoff of r c = 10 A, an Ewald parameter of 
(3 = 0.17A -1 and fifth-order interpolation (P = 5) was used to calculate the electrostatic 
potential. The grid dimension was set to 20 x 20 x 45. 

The parameters of the PPPM for dispersion are the real space cutoff, the Ewald param- 
eter, the interpolation order, and the grid spacing in each dimension. The influence of the 



different parameters is already described at the beginning of Section III Instead of exploring 
this six-dimensional parameter space, we set the interpolation order to P = 5 and the real 
space cutoff r c = 10.0 A for the hexane system. This choice of parameters was made because 
these values are commonly used in MD simulations, although they are in principle arbitrary. 
We do not claim that these are the optimal choices. For example, using the long-range dis- 
persion solver allows experimenting with smaller values for the real space cutoff and might 
in this way improve the performance of the calculations. Furthermore, the grid spacing was 
equal in all three dimensions in the simulations described below, as near cubic grids usually 
provide most accuracte calculations. 

This reduction of the parameter space allows for determining suitable simulation pa- 
rameters with less effort, but permits reaching a wide range of accuracy in either real or 
reciprocal space. As the real space cutoff is fixed, the real space accuracy depends only on 
the Ewald parameter, which is therefore used in the following simulations to tune the real 
space accuracy. In principle, we could also have fixed the Ewald parameter beforehand and 
modified the real-space cutoff in our simultations to tune the real-space accuracy, but we 
decided against it to have better control over the real space calculation time. For a given 
Ewald parameter and the other parameters fixed, the grid spacing can be altered to tune 
the reciprocal space accuracy, even though the reachable reciprocal-space accuracy is not 
unlimited for a fixed Ewald parameter. 

We performed surface tension calculations with different settings for the two remaining 
parameters, the Ewald parameter and the uniform grid spacing. We examined the resulting 
surface tensions and liquid densities. In addition we determined the rms error in the total 
forces as well as the real and reciprocal space contributions to the error by comparing 
the forces calculated for a single snapshot of an equilibrated systems to forces that were 
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FIG. 5. Surface tension and density of hexane as a function of the total error in the calculated 
forces. The arrows point in the direction of increasing Ewald parameter. 
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FIG. 6. Surface tension, density, and errors in the forces in simulations of hexane. In the lower 
left graph, the red squares correspond to the real space error, while the triangles and circles 
correspond to the reciprocal space error when using the fine grid and coarse grids, respectively. 
The circles in the lower right graph correspond to the reciprocal space error. The Ewald parameter 
is (5 = 0.28 A -1 in all figures on the right side. Dotted lines are error estimates calculated with 
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calculated using a large real-space cutoff and a very small grid spacing. 

The results for the surface tension and density of hexane are given as a function of the 
total rms error in the forces in Figure [5j In simulations with fewer grid points, the total 
error is always dominated by the reciprocal space error. In simulations with smaller grid 
spacings, where the number of grid points in the x-direction n x = 24, the real and reciprocal 
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space error are approximately equal for the highest achieved total accuracy at (3 = 0.28 A -1 . 
As can be seen from Figure [6j the real space error dominates for smaller values of the Ewald 
parameter 0, whereas the reciprocal space error dominates for larger values of (3. 

As the total error decreases, the simulated surface tensions and densities plateau, indicat- 
ing that further increases in accuracy, which can be obtained by using even finer grids and 
larger values for the Ewald parameter, will offer little benefit in the accuracy of the mea- 
sured quantities. Decreasing the Ewald parameter, thereby increasing the real space error, 
strongly influences the simulated quantities. In contrast, increasing the Ewald parameter 
and in this way increasing the reciprocal space error has less influence on the results. Phys- 
ical data begin to change for reciprocal space errors above approximately 0.01 kcal mol -1 
A -1 . For the examined quantities, the real space error has a stronger influence on the results 
than the reciprocal space error. The reason for this observation is that an increasing real 
space error leads to increasing underprediction of the cohesion of a simulated system. For 
simulations of quantities in which the cohesion does not influence the reults, the influence 
of the real and reciprocal space error will possibly be different. 

The data given in Figure [5] are also given on the left side of Figure [6] as a function of 
the Ewald parameter. These results, in combination with those from Figure [5j show that 
an Ewald parameter of approximately = 0.28 A -1 in combination with a real space cutoff 
r c = 10 A provides a sufficient real space accuracy for the performed simulations. 

As the results from Figure [5] indicate that increasing the reciprocal space error does not 
alter the obtained physical data strongly, we have performed further simulations with fixed 
Ewald parameter with varying grid spacings. Results of these simulations are given on the 
right side of Figure [6} Increasing the number of grid points n x in the x-direction beyond 
12 does not alter either the simulated density or surface tensions, although the error in the 
forces continues to decrease. However, the extended running times required for the finer 
meshes make these higher fidelity calculations computationally undesirable. 

Therefore, we choose j3 = 0.28 A -1 and the grid spacing h « 4. 17 A as these parameters 
provide sufficient accuracy. Examining the influence of the parameters of the LJ system 
provided similar results as those described above. The parameters we obtain for the LJ 
system are /3 = 1.1a -1 and h ~ 1.22cx for Interpolation order P = 5 and a real space cutoff 
of r c = 3o\ The corresponding simulations and results are described in the supporting 
information®!. 
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VII. APPLICATION OF THE SOLVER 



To compare our algorithms to existing implementations — a plain cutoff or the Ewald 
suirpH — we j-^yg performed simulations with systems of LJ particles, SPC/E water and 
hexane modeled with the OPLS all-atom force field.^ These systems cover a model system 
as well as realistic systems in which Coulomb interactions (water) and dispersion interac- 
tions (hexane) dominate. Furthermore, these systems have already been studied and allow 
comparison to results from the literature] 10 * 11 * 42 * ^ A comparison with results from Shi et 
alP^ is of special interest, as they have also used a PPPM dispersion method to determine 
the surface tension of SPC/E water. 



A. Lennard-Jones particles 

The Lennard-Jones simulations were performed in a box with volume 11.01 x 11.01 x 
176.16a 3 and 4000 particles that were placed randomly in a subvolume at the center of 
the box. After minimization using a soft potential, the system was equilibrated for 100 000 
timesteps. The timestep was set to 0.005 r, where r = a^/rafe. Simulations were exe- 
cuted at reduced temperatures T* = k B T/e G {0.7,0.85,1.1,1.2} using a Nose-Hoover^ 
thermostat with damping factor 10 r. The equations of motion were solved using a velocity 
Verlet algorithm.^ Afterwards, simulations were run for another 1 000 000 timesteps with 
the same conditions. During that time, instantaneous surface tensions were calculated every 
timestep. Configurations were stored every 1 000 timesteps to calculate the density profile. 
For simulations without a long-range dispersion solver, we examined cutoffs of 2.5a, 5cr, and 
7.5cr. Simulations with an Ewald solver were performed with cutoffs of 3a, 4a, and 5a. We 
relied on automatic generation of the Ewald parameter and the cutoff for the k vectors. We 
used the value of 0.05 as the desired relative accuracy in the forces.^ The resulting Ewald 
parameters were 0.60a -1 , 0.45a -1 , and 0.36a -1 ; the number of k vectors were 1616, 677, 
and 320 for the different cutoffs. In simulations with the dispersion PPPM we used cutoffs 
of 3a, 4a, and 5a. We used P = 5, /3 — 1.1a -1 and a grid with 9 x 9 x 144 mesh points in 



agreement with our results from Section VI 



Results are given in Table [TTJ Overall, we find good agreement with results from the 
literature] 10 * 43 * The simulated densities and surface tensions show a strong dependence on 
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TABLE II. Results of the validation runs for the LJ particles. Uncertainties given in parentheses 









Surface tension, ea 2 




T* 


solver 


r c (<t) Pliq 


[a— 3) 7 P 


It 7 


t c (ms) 


0.7 


cutoff 


2.5 


0.7865 0.588(30) 


0.327 0.915(30) 


3.7 






5.0 


0.8349 1.006(30) 


0.125 1.131(30) 


24.8 






7.5 


0.8390 1.112(30) 


0.057 1.169(30) 


80.2 




Ewald 


3.0 


0.8332 1.085(30) 


- 1.085(30) 


123.2 






4.0 


0.8371 1.121(30) 


- 1.121(30) 


77.5 






5.0 


0.8393 1.134(30) 


- 1.134(30) 


87.63 




PPPM 


3.0 


0.8404 1.158(30) 


- 1.158(30) 


19.59 






4.0 


0.8407 1.167(30) 


- 1.167(30) 


32.83 






5.0 


0.8408 1.157(30) 


- 1.157(30) 


52.67 


0.85 cutoff 


2.5 


0.6996 0.341(22) 


0.221 0.562(22) 


3.1 






5.0 


0.7672 0.700(26) 


0.098 0.798(26) 


22.6 






7.5 


0.7730 0.781(32) 


0.046 0.827(32) 


73.2 




Ewald 


3.0 


0.7651 0.742(24) 


- 0.742(24) 


122.2 






4.0 


0.7706 0.799(26) 


- 0.799(26) 


75.0 






5.0 


0.7732 0.803(28) 


- 0.803(28) 


76.0 




PPPM 


3.0 


0.7748 0.817(26) 


- 0.817(26) 


18.81 






4.0 


0.7758 0.829(24) 


- 0.829(24) 


30.00 






5.0 


0.7756 0.829(28) 


- 0.829(28) 


48.24 


1.1 


cutoff 


2.5 


n.a. 0.023(26) 


n.a. 0.023(26) 


1.4 






5.0 


0.6282 0.278(26) 


0.042 0.320(26) 


15.5 






7.5 


0.6385 0.293(24) 


0.026 0.319(24) 


50.7 




Ewald 


3.0 


0.6243 0.270(26) 


- 0.270(26) 


118.2 






4.0 


0.6354 0.293(24) 


- 0.293(24) 


68.4 






5.0 


0.6452 0.315(26) 


- 0.315(26) 


56.2 




PPPM 


3.0 


0.6451 0.314(24) 


- 0.314(24) 


15.69 






4.0 


0.6448 0.330(26) 


- 0.330(26) 


23.57 






5.0 


0.6462 0.302(22) 


- 0.302(22) 


37.19 


1.2 


cutoff 


2.5 


n.a. 0.001(20) 


n.a. 0.001(20) 


1.4 






5.0 


0.5613 0.113(20) 


0.025 0.138(20) 


11.8 






7.5 


0.5725 0.159(26) 


0.013 0.172(26) 


38.8 




Ewald 


3.0 


0.5497 0.128(24) 


- 0.128(24) 


118.9 






4.0 


0.5647 0.141(26) 


- 0.141(26) 


65.4 






5.0 


0.5728 0.159(24) 


- 0.159(24) 


49.9 




PPPM 


3.0 


0.5767 0.164(24) 


- 0.164(24) 


13.73 






4.0 


0.5757 0.155(22) 


- 0.155(22) 


21.10 






5.0 


0.5766 0.154(26) 


- 0.154(26) 


32.24 
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FIG. 7. Measured density profiles for simulations of LJ particles: (top) simulations with cutoff 
r c = 2.5a and (bottom) simulations at a reduced temperature of T* = 1.1. 

the chosen cutoff in simulation without a long-range dispersion solver. For simulations at 
higher temperatures, systems with small cutoffs were so close to the critical point that error 
functions were no longer appropriate for describing the density profile, as can be seen from 
Figure [7j Agreement between simulated data with and without long-range dispersion solver 
can only be obtained when using a large cutoff in simulations without the long-range solver. 

Unlike the simulations with a long-range cutoff, the results for the dispersion PPPM 
method do not show a dependence on the switching radius. For the Ewald sum, a slight 
dependence of the physical data remains, which we attribute to the automated parameter 
generation routine in combination with the specified accuracy. 



B. SPC/E water 

Simulations with SPC/E water were performed with 5 000 water molecules in a box of 
50 x 50 x 150 A 3 . The initial configurations of the particles were created using PackmolP^ 
If not explicitly given in the following, the simulation settings were as those for hexane 
described in Section [VTl 

Simulations were executed at 300, 350, and 400 K. For each solver, we used cutoffs of 10, 
12, and 16 A for the sum in real space for dispersion and Coulomb interaction. The SHAKE 
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algorithm^ was U sed to constrain the bond lengths and bond angles. 

A PPPM 21 solver was used for long-range electrostatics in simulations with a plain cutoff 
for dispersion interaction. We picked interpolation order P = 5 and a grid of 24 x 24 x 54 
mesh points as grid parameters. The Ewald parameter was = 0.255, 0.226, and 0.184 A -1 
for the three different cutoffs. In simulations with the traditional Ewald sum for dispersion 
and Coulomb interactions, the Ewald parameter and number of k vectors were generated 
for a desired relative accuracy of 0.05. The Ewald parameters were set to approximately 
0.18, 0.15, and 0.11 A -1 and the number of k vectors were 748, 436, and 183 for the different 
cutoffs. In simulations with a PPPM solver for dispersion, we used interpolation order 



P = 5. Following our results from Section VI, a grid with 12 x 12 x 36 mesh points was used 



for the dispersion interactions. The Ewald parameter for dispersion was set to (3 = 0.28 A -1 
for all cutoffs. The parameters used for the long-range solver for the Coulomb interaction 
were the same as those in simulations with a plain cutoff for dispersion. 

Table |III| shows the results of the simulations. When not using a long range solver, the 
simulated density shows slight dependence on the chosen cutoff radius, whereas practically 
no dependence can be observed when using a long-range solver for dispersion. For simulated 
surface tensions, neither the cutoff nor the chosen dispersion solver have a strong influence. 
The weak or non-existent influence on physical properties of the way dispersion interactions 
are calculated is due to the fact that Coulomb, and not dispersion, interactions are the 
dominant contribution to the interactions in this system. 

Again, our results are in good agreement with the majority of the literature) 10 * 42 ^ how- 
ever, they differ substantially from those reported by Shi et alP^, who performed simulations 
of SPC/E with a PPPM for dispersion, too. For example, their result for the surface tension 
at 300 K is more than 70mN/m (read from Fig. 6 in Ref. 1251) . whereas the surface tensions 
in our simulations are always about 60mN/m, consistent with other studies. To ensure 
the validity of our results we have run an additional simulation with increased accuracy, in 
which we set the Ewald parameter for dispersion to /3 = 0.3 A -1 , the interpolation order 
to P = 5 and the grid spacing to h ~ 1.56 A corresponding to 32 x 32 x 96 mesh points. 



Results of this simulation, marked with a dagger in Table III, are in good agreement with 
the rest of our results. The increased value for the surface tension in simulations by Shi et 
al. might be related to the small number of water molecules (800) in their simulation or the 
choice of the Ewald parameter (0.9, units not given), but is most likely caused by their short 
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TABLE III. Results of the validation 
was run with higher precision. 



runs for the SPC/E water. Simulation marked with a dagger 



Surface tension, mN/m 
T (K) solver r c (A) p Uq (kg/L) -y p -y t 7 t c (ms) 



300 


cutoff 


10.0 


0.9882 54.59(100) 


5.27 59.86(100) 


169 






12.0 


0.9918 56.38(100) 


3.72 


60.10(100) 


259 






16.0 


0.9944 57.51(84) 


2.12 


59.63(84) 


531 




Ewald 


10.0 


0.9964 59.39(100) 




59.39(100) 


455 






12.0 


0.9962 60.84(100) 


_ 


60.84(100) 


474 






16.0 


0.9967 60.56(100) 


_ 


60.56(100) 


756 




PPPM 


10.0 


0.9965 60.72(90) 


_ 


60.72(90) 


229 






12.0 


0.9964 60.11(80) 


- 


60.11(80) 


364 






16.0 


0.9963 59.64(90) 


_ 


59.64(90) 


737 






10.0 


0.9964 61.06(80) 


_ 


61.06(80) f 


- 


350 


cutoff 


10.0 


0.9539 47.40(60) 


4.81 


52.21(60) 


166 






12.0 


0.9576 48.71(60) 


3.42 


52.13(60) 


254 






16.0 


0.9607 49.78(70) 


1.96 


51.74(70) 


517 




Ewald 


10.0 


0.9617 52.55(80) 


_ 


52.55(80) 


448 






12.0 


0.9622 51.72(70) 


_ 


51.72(70) 


531 






16.0 


0.9626 52.13(70) 




52.13(70) 


743 




PPPM 


10.0 


0.9629 53.29(60) 


_ 


53.29(60) 


239 






12.0 


0.9631 52.35(70) 


_ 


52.35(70) 


350 






16.0 


0.9630 52.30(70) 


_ 


52.30(70) 


737 


400 


cutoff 


10.0 


0.9067 39.89(60) 


4.20 


44.09(60) 


165 






12.0 


0.9114 40.42(60) 


3.02 


43.44(60) 


244 






16.0 


0.9151 41.36(58) 


1.75 


43.11(58) 


494 




Ewald 


10.0 


0.9164 42.97(60) 




42.97(60) 


438 






12.0 


0.9168 43.71(70) 




43.71(70) 


512 






16.0 


0.9172 43.34(60) 




43.34(60) 


716 




PPPM 


10.0 


0.9177 43.98(60) 




43.98(60) 


240 






12.0 


0.9178 43.89(60) 




43.89(60) 


345 






16.0 


0.9178 43.44(60) 




43.44(60) 


710 
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sampling time of only 100 000 timesteps, as substantially longer run times are required to 
achieve equilibration for water at an interface.^ 



C. Hexane 



If not given in the following, all settings for the hexane simulations were as those reported 



in Section [VI} We studied temperatures of 300, 350, and 400 K and cutoffs of 10, 12, and 
16 A. 

In simulations with a plain cutoff for dispersion, a PPPM solver was used for electrostat- 
ics. The grid dimension was set to 20 x 20 x 45 and the interpolation order to P = 5. The 
Ewald parameter was approximately /3 = 0.17, 0.16, and 0.14 A" 1 for the different cutoffs. 
The desired precision was set to 0.05 in simulations with the Ewald method for dispersion 
and Coulomb interactions. The resulting Ewald parameters and number of k vectors were 
the same as those in simulations with SPC/E water. In simulations with the PPPM for 
dispersion, the interpolation order was set to P = 5, the grid size was set to 12 x 12 x 36 
and the Ewald parameter was set to /3 — 0.28 A^ 1 in all simulations. Coulomb interactions 
were treated in a same way as in simulations with a cutoff for dispersion. 



The results are summarized in Table IV The chosen cutoff radius has a strong influence 
on the results in simulations without a long-range dispersion solver. In contrast, the results 
for Ewald summation show weak and the results for the PPPM show no dependence on 
the chosen cutoff radius. Our results for the PPPM are in good agreement with those from 
Ismail et al.^ in simulations with an Ewald sum for dispersion. This, and the fact that the 
chosen cutoff does not influence the results, confirms the validity of our simulations and the 
good choice of the Ewald and grid parameters. Our simulations with Ewald sums provide 
lower surface tensions, which is caused by insufficient accuracy in these simulations. As can 
be seen from Figure |8j the simulation results when not using a long-range dispersion solver 
approach those obtained with the PPPM when increasing the cutoff. However, even those 
with a cutoff of 16 A provide surface tensions and densities that are below those obtained 
from PPPM simulations. 
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FIG. 8. Surface tensions and densities simulated when using a PPPM or a plain cutoff for disper- 
sion. Results obtained when not using a long-range solver strongly depend on the chosen cutoff 
and approach the results of the simulation with the PPPM with increasing cutoff size. 

VIII. PERFORMANCE COMPARISON 



To measure the simulation time, each of the simulations in Section VII was run for 1 000 



timesteps on a single core. The resulting computation times per timestep t c are given in the 
last column of Tables [IT] to IV The simulations were executed on Intel Harpertown E5454 
processors with eight 3.0 GHz Xeon cores. 

For a fair comparison between the different solvers, one should consider different solvers 
at the same temperature with the cutoff that provides the results that are obtained in the 
limit of high accuracy simulations. If for a given solver and temperatures different cutoffs 
provided accurate results, the fastest of those simulations should be used. 

A quick comparison of the PPPM with the Ewald shows that simulations with the PPPM 
were faster in all cases. As the Ewald sum was always slower than the PPPM, we omit 
comparisons between the Ewald solver and the plain cutoff and continue with comparing 
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TABLE IV. Results of the validation runs for the OPLS hexane 



Surface tension, mN/m 
T (K) solver r c (A) p liq (kg/L) -y p j t 7 t c (ms) 



300 


cutoff 


10.0 


0.6058 7.83(50) 


4.59 12.42(50) 


137 






12.0 


0.6251 10.21(50) 


3.75 13.96(50) 


207 






16.0 


0.6367 13.00(50) 


2.34 15.34(50) 


421 




Ewald 


10.0 


0.6368 14.40(50) 


- 14.40(50) 


409 






12.0 


0.6385 14.91(50) 


- 14.91(50) 


414 






16.0 


0.6410 14.91(56) 


- 14.91(56) 


634 




PPPM 


10.0 


0.6434 16.41(50) 


- 16.41(50) 


201 






12.0 


0.6439 16.16(50) 


- 16.16(50) 


352 






16.0 


0.6453 15.89(50) 


- 15.89(50) 


601 


350 


cutoff 


10.0 


0.5237 2.18(40) 


2.13 4.31(40) 


118 






12.0 


0.5534 4.78(40) 


2.20 6.98(40) 


182 






16.0 


0.5721 7.44(50) 


1.71 9.15(50) 


383 




Ewald 


10.0 


0.5722 8.09(50) 


- 8.09(50) 


384 






12.0 


0.5778 8.52(40) 


- 8.52(40) 


438 






16.0 


0.5805 9.03(40) 


- 9.03(40) 


575 




PPPM 


10.0 


0.5823 9.97(44) 


- 9.97(44) 


183 






12.0 


0.5839 9.77(60) 


- 9.77(60) 


276 






16.0 


0.5851 9.89(44) 


- 9.89(44) 


547 


400 


cutoff 


10.0 


n.a -1.55(32) 


n.a. -1.55(32) 


92 






12.0 


0.4467 0.30(36) 


0.74 1.04(36) 


149 






16.0 


0.4905 1.83(36) 


0.85 2.68(36) 


316 




Ewald 


10.0 


0.4881 2.26(50) 


- 2.26(50) 


368 






12.0 


0.5037 3.07(32) 


- 3.07(16) 


402 






16.0 


0.5039 3.43(40) 


- 3.43(40) 


498 




PPPM 


10.0 


0.5099 4.59(36) 


- 4.59(36) 


158 






12.0 


0.5106 4.66(40) 


- 4.66(40) 


234 






16.0 


0.5157 4.45(46) 


- 4.45(46) 


452 
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simulations with a PPPM to those with a plain cutoff. For the LJ system, accurate surface 
tensions and densities were obtained only in simulation in which the cutoff was r c = 7.5a 
in simulations without long-range dispersion solvers. For simulations with the PPPM for 
dispersion r c = 3<r provides proper results with maximum efficiency Comparing the cor- 
responding simulation times shows the computational superiority of using the PPPM for 
dispersion. 

In simulations with water, the results of the comparison are different when examining the 
simulated surface tensions or simulated densities. As surface tensions were approximately 
the same throughout all simulations, the PPPM was outperformed by the plain cutoff for 
this comparison. The reason for the observed behavior is not the choice of parameters, 
but the dominance of Coulomb interaction in the examined system. When comparing the 
simulation times it has to be kept in mind that the simulations with a plain cutoff are 
incorrect during the simulations and are only corrected a posteriori. If this correction is not 
possible after the simulations, or if a correct value of the surface tension is required during 
the simulations for any reason, a larger cutoff is required in simulations with a plain cutoff. 
Simulations with a PPPM are more efficient in such cases. The increase in simulation time 
is about 35 percent when using the PPPM for dispersion compared to a plain cutoff. 

If highly accurate simulated densities are important, then the cutoff for dispersion interac- 
tions should be chosen to be at least r c = 12 A in simulations without long-range dispersion 
solvers. Comparing the simulation time of these simulations to those with a PPPM with a 
cutoff r c = 10 A shows that using the PPPM is computationally favorable in this case. 

For simulations of hexane, in which dispersion interactions dominate, the largest cutoff 
has to be selected in simulations without long-range dispersion solvers, whereas a small 
cutoff is sufficient in simulations with the PPPM. As a consequence, the simulations with 
the PPPM were much faster than those without a long-range dispersion solver. 

We would like to note that simulations with the long-range dispersion solver were run 
without tabulating the pair potential. Including this tabulation will results in additional 
speed-up of the simulations and might change the comparison above. 
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IX. CONCLUSIONS AND OUTLOOK 



We present a PPPM algorithm for dispersion interactions that calculates long-range dis- 
persion forces accurately and efficiently for inhomogeneous systems. When used correctly, 
this solver strongly reduces the error that is caused when truncating the pair potential at a 
plain cutoff and thus provides a better description of the physics of a simulated system. 

We derived and tested error estimates that describe the effect of the parameters of the 
PPPM on simulation results. The presented error estimates are only valid for bulk phase 
systems and should not be relied on when simulated energies or pressures are relevant. 

For not having to rely on the presented error measures, we explored the parameter space 
to provide parameters that can be used in future simulations. For real physical systems 
of surfaces, a combination of the Ewald parameter (5 = 0.28 A -1 , the interpolation order 
P = 5, the grid spacing h ~ 4.17 A, and the real space cutoff r c = 10.0 A provided good 
results for different materials at different temperatures. 

We have applied the developed algorithm to study the surface tension of LJ particles, 
SPC/E water, and hexane. The described algorithm outperforms Ewald summation for all 
systems that were examined in this study in terms of simulation time. 

Comparing the PPPM to simulations with a plain cutoff show that the PPPM is favorable 
when correction methods cannot be applied or correction methods do not work properly, as 
for example near the critical point or in some of our hexane and LJ simulations. In systems 
that are dominated by dispersion, the PPPM outperforms simulations without long-range 
solvers in terns of computation time, as latter simulations require larger cutoffs. 

For strongly charged systems, the PPPM provides a benefit in simulation time only if 
densities are needed at a high accuracy. In other cases a plain cutoff is favorable in terms 
of computation time here. However, the advantage of correctness during the simulation and 
the applicability to arbitrarily shaped surfaces remains. 
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